#include <iostream>
// #include <cstdio>
#include <stdio.h>
#include <string>

#ifdef _WIN32
#include <direct.h>

#elif __APPLE__ || __linux__
#include<unistd.h>
// 使用armadillo矩阵库
#include <armadillo>
#endif

#include <memory>
#include <cmath>

#include "EmoRecBase.h"
#include "L1WeightMatrix.h"
#include "ColorfulLog.h"
#include "EmotionFilter.h"
#include "TrainMethod.h"

#define GET_TIME


// #define DEBUG
// 构造函数
// 只用了63通道，训练时去掉EOG眼电通道
/**********************************
 * 函数名：EmoRecBase
 * 功能：构造函数
 * 输入：无
 * 输出：无
 * 备注：无
**********************************/
EmoRecBase::EmoRecBase(ConfigImport* sysconf):
    sys_config(*sysconf), network_cache(sys_config.getMapVal("Cache_path"),1,0)
{
    trainMethod_ptr = new SVM(sys_config.getMapVal("SVM_modelPath"));
    // 初始化缓存策略
    cacheConfig_init();
}

EmoRecBase::EmoRecBase(ConfigImport &configPara)
    :sys_config(configPara),network_cache(configPara.getMapVal("Cache_path"),1,0)
{
    trainMethod_ptr = new SVM(sys_config.getMapVal("SVM_modelPath"));
}

EmoRecBase::~EmoRecBase(){
    if(trainMethod_ptr != nullptr){
        delete trainMethod_ptr;
        trainMethod_ptr = nullptr;
        std::cout << "In EmoRecBase: trainMethod_ptr is deleted .." << std::endl;
    }
}

/**********************************
 * 函数名：generateLabelFromFile
 * 功能：从文件中导入标签，用于离线数据分析
 * 输入：fileName: 文件名
 *      labelType: 0:train_label, 1:test_label
 * 输出：true：成功，false：失败
 * 备注：无
**********************************/
bool EmoRecBase::generateLabelFromFile(std::string fileName, int labelType){
    std::ifstream labelFile_fs(fileName,std::ios::in);
    if(!labelFile_fs.is_open()){
        LOG_ERR("file \"%s\" is not exit, open failed.\n",fileName.c_str());
        return false;
    }
    // // 读取文件大小
    labelFile_fs.seekg(0, std::ios::end);
    std::streampos fileSize = labelFile_fs.tellg();
    labelFile_fs.seekg(0, std::ios::beg);
    // 计算float数组的元素数量
    size_t floatCount = fileSize / sizeof(float);
    // 创建存储float数据的数组
    std::vector<float> floatArray(floatCount);
    // 从文件中读取数据到数组中
    labelFile_fs.read(reinterpret_cast<char*>(floatArray.data()), floatCount*sizeof(float));
    // 关闭文件
    labelFile_fs.close();
    // 根据数据生成标签
    switch (labelType)
    {
        case 0:
            // train_label
            for(int i=0; i<floatArray.size(); ++i){
                if(train_label.empty()){
                    train_label = std::vector<int>(floatArray[i],(i/3)+1);
                    continue;
                }
                std::vector<int> label_temp = std::vector<int>(floatArray[i],(i/3)+1);
                train_label.insert(train_label.end(),label_temp.begin(),label_temp.end());
            }
            break;
            
        case 1:
            // test_label
            for(int i=0; i<floatArray.size(); ++i){
            if(test_label.empty()){
                test_label = std::vector<int>(floatArray[i],i+1);
                continue;
            }
            std::vector<int> label_temp = std::vector<int>(floatArray[i],i+1);
            test_label.insert(test_label.end(),label_temp.begin(),label_temp.end());
        }
            break;
        
        default:
            break;
    }
    return true;
}

/**********************************
 * 函数名：cacheConfig_init
 * 功能：
 * 输入：无
 * 输出：无
 * 备注：无
**********************************/
void EmoRecBase::cacheConfig_init(){
    // 全局不使用缓存功能
    if(!sys_config.get_GL_Cache()){
        networkCache_prepared = false;
        refresh_networkCache = false;
        return;
    }
    // 刷新所有缓存
    if(sys_config.get_CL_Cache()){
        networkCache_prepared = false;
        refresh_networkCache = true;
        return;
    }

    // 刷新脑网络缓存
    if(sys_config.is_cor_cache_refresh()){
        networkCache_prepared = false;
        refresh_networkCache = true;
        return;
    }

    // 使能脑网络缓存且不刷新脑网络缓存
    if(sys_config.is_cor_cache_enable()){
        uint8_t res = network_cache.load(Train_Feature);
        if(!res){
            networkCache_prepared = false;
            refresh_networkCache = true;
            return;
        }
        networkCache_prepared = true;
        refresh_networkCache = false;
        return;
    }

    // 不使用脑网络缓存
    networkCache_prepared = false;
    refresh_networkCache = false;
}

/**********************************
 * 函数名：train
 * 功能：训练模型
 * 输入：无
 * 输出：0:训练成功, 
 *      1:网络矩阵为空
 *      2:标签为空
 *      3:标签长度和数据长度不对齐
 * 备注：无
**********************************/
int EmoRecBase::train(){
    printf("begin trainning...\n");
    if(Train_Feature.n_elem == 0){
        printf("Please run \"network_created(arma::mat &data, arma::mat &network_output)\" to creat a network.\n");   
        return 1;
    }
    if(!train_label.size()){
        printf("train label is empty!\n");
        return 2;
    }
    arma::vec train_label_vec = arma::conv_to<arma::vec>::from(train_label);
    if(train_label_vec.n_elem != Train_Feature.n_rows){
        LOG_ERR("train_label_vec.n_elem = %lld,Train_Feature.n_rows = %lld\n", train_label_vec.n_elem, Train_Feature.n_rows);
        return 3;
    }

    #ifdef GET_TIME
    auto BT_LSL1 = std::chrono::steady_clock::now();
    #endif
    LaplaceSupervisedL1(Train_Feature, train_label_vec, l1NormEigvector);
    Train_Feature = Train_Feature * l1NormEigvector;
    #ifdef DEBUG
    std::cout << Train_Feature.n_rows << "x" << Train_Feature.n_cols << std::endl;
    // Train_Feature.print("Train_Feature:");
    #endif // DEBUG
    
    #ifdef GET_TIME
    auto AT_LSL1 = std::chrono::steady_clock::now();
    double duration_millsecond = std::chrono::duration<double, std::milli>(AT_LSL1 - BT_LSL1).count();
	std::cout << "LSL1_times :" << duration_millsecond << "ms" << std::endl; 
    #endif
    // 创建一个分类模型
    trainMethod_ptr->create_model(Train_Feature, train_label_vec);

    // 保存特征向量
    l1NormEigvector.save("../Cache/normEigvector.cache",arma::raw_ascii);
    return 0;
}

/**********************************
 * 函数名：predict
 * 功能：基于现有的模型，预测当前一个脑网络的标签
 * 输入：predict_data：待预测的矩阵数据
 * 输出：无
 * 备注：无
**********************************/
int EmoRecBase::predict(arma::mat &predict_data){
    arma::mat pre_network;
    // network_created(predict_data, pre_network);
    // // pre_network.print("pre_net:");
    // if(pre_network.n_elem == 0){
	//     std::cerr << "pre_network is empty!!! " << std::endl;
    //     return 3;
    // }
    if(l1NormEigvector.n_elem == 0){
        l1NormEigvector.load("../Cache/normEigvector.cache",arma::raw_ascii);
    }
    pre_network = predict_data;
    pre_network *= l1NormEigvector;
	// std::cout << "pre_network:" << pre_network.n_rows << "x" << pre_network.n_cols << std::endl;
    return trainMethod_ptr->data_predict(pre_network);
}

/**********************************
 * 函数名：insertLabel
 * 功能：追加训练标签
 * 输入：无
 * 输出：当前标签长度
 * 备注：无
**********************************/
int EmoRecBase::insertLabel(int inputLabel){
    train_label.push_back(inputLabel);
    return train_label.size();
}

int EmoRecBase::insertLabel(std::vector<int>& inputLabel){
    // 这里可能会拖慢解码时间
    train_label.insert(train_label.end(),inputLabel.begin(),inputLabel.end());
    return train_label.size();
}


void EmoRecBase::network_created(arma::cube &inputData){
    if(!inputData.n_elem){
        std::cerr << "In EmoRecBase::network_created,inputData is empty!" << std::endl;
        return;
    }
    // 创建滤波器
    EmotionFilter emotionfilter = EmotionFilter(sys_config.getMapVal("EmotionFilterPara"));
    int times = inputData.n_slices;
    for(int trial_index=0; trial_index<times; ++trial_index){
        #ifdef GET_TIME
        auto BT_filt = std::chrono::steady_clock::now();
        #endif // GET_TIME
        
        // 取出一个trial
        arma::mat n_trial;
        n_trial = inputData.slice(trial_index);
        // 滤波
        arma::mat outMat(n_trial.n_rows,n_trial.n_cols);
        emotionfilter.filt(n_trial,outMat);
        #ifdef GET_TIME
        auto AT_filt = std::chrono::steady_clock::now();
        #endif // GET_TIME

        // 第二步，相关网络计算
        arma::mat cor_network;
        #ifdef GET_TIME
        auto BT_cor = std::chrono::steady_clock::now();
        #endif
        creat_network_cor(outMat, cor_network);
        #ifdef GET_TIME
        auto AT_cor = std::chrono::steady_clock::now();
        #endif

        // 将脑网络拉伸成一个行向量
        arma::mat unitUpperTri = arma::ones<arma::mat>(n_trial.n_rows, n_trial.n_rows);
        // 将下三角部分置为 0
        unitUpperTri = arma::trimatu(unitUpperTri);
        unitUpperTri.diag().zeros();
        arma::mat resMat = cor_network(arma::find(unitUpperTri));
        arma::rowvec rowVec = arma::vectorise(resMat).t();

        // // 去心
        // cor_network.diag().zeros();
        // // 提取上三角矩阵
        // arma::mat upperTri = trimatu(cor_network);
        // // upperTri.print("upperTri:");
        // // bool res = upperTri.save("../Data/upperTri.txt",arma::arma_ascii);
        // // // 直接将上三角矩阵的非零元素存储到行向量中
        // arma::rowvec rowVec = arma::nonzeros(upperTri).t();
        
        static int once = 1;
        if(once){
            once = 0;
        }
        else{
            if(rowVec.n_cols != Train_Feature.n_cols){
                printf("rowVec.n_cols = %lld,Train_Feature.n_cols = %lld\n", rowVec.n_cols, Train_Feature.n_cols);
                printf("trial_index/times : %d/%d\n",trial_index,times);
                printf("Train_Feature.n_rows=%lld\n",Train_Feature.n_rows);
                cor_network.print("cor_network:");
                rowVec.print("rowVec:");
                while(1);
            }
        }
        Train_Feature = arma::join_vert(Train_Feature,rowVec);
    }
    // printf("cubeLen:%lld\n",inputData.n_slices);
    // Train_Feature.print("Train_Feature:");
}



// =============================================================================== //
void EmoRecBase::train(arma::cube &data_input){
    data_org = data_input;
    // generated_label(sys_config.getMapVal("labelFileName"));
    generateLabelFromFile(sys_config.getMapVal("labelFileName"),0);
    // 如果有缓存，则不需要重建网络
    if(sys_config.get_GL_Cache()){
        if(!sys_config.get_CL_Cache()){
            if(sys_config.is_cor_cache_enable()){
                if(!sys_config.is_cor_cache_refresh()){
                    uint8_t res = network_cache.load(Train_Feature);
                    if(!res){
                        network_created();
                        network_cache.save(Train_Feature);
                        #ifdef DEBUG
                        std::cout << "network saving ..." << std::endl;
                        #endif
                    }

                }
                else{
                    network_created();
                    network_cache.save(Train_Feature);
                    #ifdef DEBUG
                    std::cout << "network saving ..." << std::endl;
                    #endif
                    
                }
            }
            else{
                network_created();
                #ifdef DEBUG
                std::cout << "network refresh ..." << std::endl;
                #endif

            }
        }
        else{
            network_created();
            #ifdef DEBUG
            std::cout << "network refresh ..." << std::endl;
            #endif
            if(sys_config.is_cor_cache_enable()){
                network_cache.save(Train_Feature);
                #ifdef DEBUG
                std::cout << "network saving ..." << std::endl;
                #endif
                
            }
        }
    }
    else{
        network_created();
        #ifdef DEBUG
        std::cout << "network refresh ..." << std::endl;
        #endif
    }

    arma::vec train_label_vec = arma::conv_to<arma::vec>::from(train_label);


    #ifdef GET_TIME
    auto BT_LSL1 = std::chrono::steady_clock::now();
    #endif
    LaplaceSupervisedL1(Train_Feature,train_label_vec,l1NormEigvector);
    Train_Feature = Train_Feature * l1NormEigvector;
    // std::cout << Train_Feature.n_rows << "x" << Train_Feature.n_cols << std::endl;
    // Train_Feature.print("Train_Feature:");

    #ifdef GET_TIME
    auto AT_LSL1 = std::chrono::steady_clock::now();
    #endif
    double duration_millsecond = std::chrono::duration<double, std::milli>(AT_LSL1 - BT_LSL1).count();
	std::cout << "LSL1_times :" << duration_millsecond << "ms" << std::endl; 
    
    // 创建一个分类模型
    trainMethod_ptr->create_model(Train_Feature, train_label_vec);
    
}

void EmoRecBase::predict_test(arma::cube &predict_data){
    // generated_test_label(sys_config.getMapVal("test_labelFileName"));
    generateLabelFromFile(sys_config.getMapVal("test_labelFileName"),1);

    // BTMethod* trainMethod_ptr = new SVM("../Data/svm_model");
    int predict_label = 0;
    int correct = 0;
    int total = 0;

    for (int i=0; i<predict_data.n_slices; ++i){
        arma::mat pre_network;
        network_created(predict_data.slice(i), pre_network);
        pre_network *= l1NormEigvector;
        predict_label = trainMethod_ptr->data_predict(pre_network);
        if(predict_label == test_label[i])
			++correct;
        ++total;
    }

    printf("Accuracy = %g%% (%d/%d) (classification)\n",
			(double)correct/total*100,correct,total);

	// std::cout << "pre_network:" << pre_network.n_rows << "x" << pre_network.n_cols << std::endl;
    // pre_network.print("pre_net_*:");
}

/**********************************
 * 函数名：network_cor
 * 功能：输入一个trial,创建一个相关网络
 * 输入：
 * 输出：无
 * 备注：无
**********************************/
void EmoRecBase::creat_network_cor(arma::mat &data_input, arma::mat &network_output){
    // cor方法计算的是矩阵的列之间的相关性,所以要先转置,这里要找出行之间的相关性
    network_output = arma::cor(data_input.t());
    // 输出成txt文件，用matlab检验结果
    // bool res = cor_mat.save("../Data/martix_arma.txt",arma::arma_ascii);
    // cor_mat.print();

}

/**********************************
 * 函数名：network_created
 * 功能：离线测试的train_data的脑网络创建
 * 输入：
 * 输出：无
 * 备注：无
**********************************/
void EmoRecBase::network_created(){
    // m_cor_network.clear();
    // 创建滤波器
    EmotionFilter emotionfilter = EmotionFilter(sys_config.getMapVal("EmotionFilterPara"));
    int times = sys_config.dim[2];
    for(int trial_index=0; trial_index<times; ++trial_index){
    // for(int trial_index=0; trial_index<10; ++trial_index){
        #ifdef GET_TIME
        auto BT_filt = std::chrono::steady_clock::now();
        #endif // GET_TIME
        
        // 取出一个trial
        arma::mat n_trial;
        // if(data_org.n_slices >= 468){
        //     n_trial = data_org.slice(trial_index);
        // }
        // else{
        //     n_trial = data_org.slice(0);
        // }

        n_trial = data_org.slice(trial_index);

        // 滤波
        arma::mat outMat(n_trial.n_rows,n_trial.n_cols);
        emotionfilter.filt(n_trial,outMat);
        // outMat.print();
        #ifdef GET_TIME
        auto AT_filt = std::chrono::steady_clock::now();
        #endif // GET_TIME

        // 第二步，相关网络计算
        arma::mat cor_network;
        #ifdef GET_TIME
        auto BT_cor = std::chrono::steady_clock::now();
        #endif
        creat_network_cor(outMat, cor_network);
        #ifdef GET_TIME
        auto AT_cor = std::chrono::steady_clock::now();
        #endif

        // 将脑网络拉伸成一个行向量
        // 去心
        cor_network.diag().zeros();
        // 提取上三角矩阵
        arma::mat upperTri = trimatu(cor_network);
        // upperTri.print("upperTri:");
        // bool res = upperTri.save("../Data/upperTri.txt",arma::arma_ascii);
        // // 直接将上三角矩阵的非零元素存储到行向量中
        arma::rowvec rowVec = arma::nonzeros(upperTri).t();
        Train_Feature = arma::join_vert(Train_Feature,rowVec);
    }
}

/**********************************
 * 函数名：network_created
 * 功能：计算单个矩阵的脑网络，测试用
 * 输入：data:数据，network_output:脑网络接收参数
 * 输出：无
 * 备注：无
**********************************/
void EmoRecBase::network_created(arma::mat &data, arma::mat &network_output){
    // 创建滤波器
    EmotionFilter emotionfilter = EmotionFilter(sys_config.getMapVal("EmotionFilterPara"));
        #ifdef GET_TIME
        auto BT_filt = std::chrono::steady_clock::now();
        #endif // GET_TIME
        
        // 滤波
        arma::mat outMat(data.n_rows,data.n_cols);
        emotionfilter.filt(data,outMat);
        // outMat.print();
        #ifdef GET_TIME
        auto AT_filt = std::chrono::steady_clock::now();
        #endif // GET_TIME

        // 第二步，相关网络计算
        arma::mat cor_network;
        #ifdef GET_TIME
        auto BT_cor = std::chrono::steady_clock::now();
        #endif
        creat_network_cor(outMat, cor_network);
        #ifdef GET_TIME
        auto AT_cor = std::chrono::steady_clock::now();
        #endif

        // 将脑网络拉伸成一个行向量
        // 去心
        cor_network.diag().zeros();
        // 提取上三角矩阵
        arma::mat upperTri = trimatu(cor_network);
        // upperTri.print("upperTri:");
        // bool res = upperTri.save("../Data/upperTri.txt",arma::arma_ascii);
        // // 直接将上三角矩阵的非零元素存储到行向量中
        network_output = arma::nonzeros(upperTri).t();    
}


/**********************************
 * 函数名：EmoRecBase::LaplaceSupervisedL1
 * 功能：计算L1Norm 特征向量
 * 输入：InputData：mat类型，每一行表示一个网络(从实对称脑网络取上三角，去掉对角线，按列展开拉伸成一个行向量)
 *                          一共有trial数那么多列
 *      InputLabel：数据标签向量，每一行对应一个标签
 *      Eigvector：特征向量
 * 输出：无
 * 备注：无
**********************************/
int EmoRecBase::LaplaceSupervisedL1(
    arma::mat &InputData, arma::vec &InputLabel, 
    arma::mat &Eigvector){
    // 相关参数设置
    int ReducedDim = 15;
    int KNum = 20;
    L1WeightMatrix_para_typeDef options;
    // options.Metric = "Euclidean";
    // options.PCARatio = 0.99;
    // options.NeighborMode = "KNN";
    // options.k = KNum;
    // options.WeightMode = "HeatKernel";
    // options.t = 3;

    // 引入标签信息，构造类内类间邻接矩阵W
    // Index 原向量元素在排序后的向量的位置(升序)
    // arma::rowvec sort_label = sort(InputLabel, "ascend");
    arma::vec SortedLabel = sort(InputLabel, "ascend");
    arma::uvec sort_label_index = sort_index(InputLabel, "ascend");

    arma::mat SortedData = InputData.rows(sort_label_index);
    // SortedData.print("sortData:");

    // 标签的种类
    arma::vec LabelKind = arma::unique(SortedLabel);
    // 标签的个数
    unsigned int SortedLabelNum = SortedLabel.n_elem;

    // 取出类内和类间的判断矩阵
    arma::mat IndexMatrix;  // 初始为空矩阵
    for (unsigned int i = 0; i < LabelKind.n_elem; ++i) {
        // 选中当前标签值
        double SelectedLabel = LabelKind(i);
        // 获取选中标签在排序后标签列表的序号
        arma::uvec SelectedIndex = arma::find(SortedLabel == SelectedLabel);
        // 当前标签的样本个数，作为序号矩阵的列数，trial数作为矩阵的行数
        unsigned int SampleNum = SelectedIndex.n_elem;
        
        // 创建一个向量，相当于matlab的SelectedIndex
        arma::vec SelectedIndex_vector(SortedLabelNum,1,arma::fill::zeros);
        SelectedIndex_vector.elem(SelectedIndex).fill(1.0);

        // 扩充SelectedIndex_vector
        arma::mat IndexMatrix1(SortedLabelNum, SampleNum, arma::fill::zeros);
        IndexMatrix1.each_col() = arma::conv_to<arma::vec>::from(SelectedIndex_vector);

        // 附加到 IndexMatrix
        if (IndexMatrix.n_elem == 0) {
            IndexMatrix = IndexMatrix1;
        } else {
            IndexMatrix = arma::join_horiz(IndexMatrix, IndexMatrix1);
        }
    }
    // IndexMatrix.print("IndexMatrix:");
    arma::umat IndexMatrix2_bool = (IndexMatrix == 0);
    arma::mat IndexMatrix2 = arma::conv_to<arma::mat>::from(IndexMatrix2_bool);
    // IndexMatrix2.print("IndexMatrix2:");
    // IndexMatrix2_bool.print("IndexMatrix2_bool:");

    // IndexMatrix2 = IndexMatrix;
    // for(auto it = IndexMatrix2.begin(); it != IndexMatrix2.end(); ++it){
    //     if(it == 0){
    //         *it = 1;
    //     }
    //     else{
    //         *it = 0;
    //     }
    // }

    // 传参时，SortedData为稠密矩阵
    // 网络数据预处理，减去均值
    arma::rowvec rvec_sampleMean = arma::mean(SortedData);
    arma::mat meanMatrix(SortedLabelNum, rvec_sampleMean.n_elem, arma::fill::zeros); 
    meanMatrix.each_row() = arma::conv_to<arma::rowvec>::from(rvec_sampleMean);
    // meanMatrix.print("meanMatrix:");
    temp_data = arma::mat(SortedData - meanMatrix);
    // SortedData.print("SortedData:");
    // data.print("data:");

    // 构建邻接矩阵
    std::shared_ptr <L1WeightMatrix> W;
    // L1WeightMatrix W;
    
    // L1WeightMatrix中参数初始化失败会抛出异常,这里使用异常捕获来创建对象
    try{
        // KNN
        // W = L1WeightMatrix(data,std::string("Euclidean"), 0.99, L1WM::Supervised, 20, L1WM::HeatKernel,3);
        W = std::make_shared<L1WeightMatrix> (temp_data,std::string("Euclidean"), 0.99, 
                                                L1WM::KNN, 20, L1WM::HeatKernel,3);

    }
    catch(const char* msg){
        std::cout << msg << std::endl;
    }

    if(!W){
        LOG_ERR("L1WeightMatrix init error, please check the parameter!\n");
        return -1;
    }
    // 在L1WeightMatrix构造时已经运算过了,这里省略
    // arma::sp_mat W1 = W->get_weightMatrix_sp();
    // W1 = arma::max(W1,W1.t());

    // 求列和后对角化

    arma::mat FullW = W->get_weightMatrix_full();
    arma::mat D = arma::diagmat(arma::sum(FullW));
    
    // 将邻接矩阵按照标签，求出类内及类间邻接矩阵
    // 元素对应相乘，取出标签对应为1的元素
    // 逐个元素相乘
    arma::mat ClassIn = (FullW % IndexMatrix);
    arma::mat ClassBe = (FullW % IndexMatrix2);

    #ifdef DEBUG
    IndexMatrix.print("IndexMatrix:");
    IndexMatrix2.print("IndexMatrix2:");
    ClassIn.print("ClassIn:");
    ClassBe.print("ClassBe:");
    #endif // DEBUG
    
    // 计算类内类间邻接矩阵的1/2次项
    arma::mat LIn = D - ClassIn;
    arma::vec S;
    arma::mat U,V;
    arma::svd(U,S,V,LIn);
    #ifdef DEBUG
    U.print("U:");
    S.print("S:");
    V.print("V:");
    LIn.print("LIn:");
    #endif // DEBUG
    arma::mat SS = arma::pow(S,0.5);
    arma::mat SqurClassIn = U*arma::diagmat(SS);

    arma::mat LBe = D - ClassBe;
    arma::svd(U,S,V,LBe);
    #ifdef DEBUG
    U.print("U:");
    S.print("S:");
    V.print("V:");
    LBe.print("LBe:");
    #endif // DEBUG

    SS = arma::pow(S,0.5);
    arma::mat SqurClassBe = U*arma::diagmat(SS);

    #ifdef DEBUG
    SqurClassIn.print("SqurClassIn:");
    SqurClassBe.print("SqurClassBe:");
    // 重构验证SVD分解准确性
    arma::mat dig_S = arma::diagmat(S);
    arma::mat back = U*dig_S*V.t();
    back.print("back:");
    #endif // DEBUG
    // 计算目标函数初始矩阵
    arma::mat Sb = temp_data.t()*SqurClassBe;
    arma::mat Sw = temp_data.t()*SqurClassIn;

    #ifdef DEBUG
    Sb.print("Sb:");
    Sw.print("Sw:");
    #endif // DEBUG
    
    //  迭代求解正交映射向量
    arma::mat IterVector;
    // 列数等于二重循环的次数+1
    arma::mat W_Record(temp_data.n_cols,101,arma::fill::zeros);
    double eta = 0.00001;
    // 迭代参数
    arma::rowvec gamma = {-10,-1,-0.1,-0.001,-0.0001,0.0001,0.001,0.01,0.1,1,10};
    arma::rowvec J = arma::zeros(1,101);

    #ifdef DEBUG
    gamma.print("gam:");
    #endif // DEBUG

    // 传出向量若为空，则进行初始化
    if(Eigvector.is_empty()){
        Eigvector = arma::mat(temp_data.n_cols,ReducedDim,arma::fill::zeros);
    }

    for(int Iter=0; Iter<ReducedDim; ++Iter){
        // 随机化定义迭代起始向量,正态分布
        arma::mat Wvector = arma::randn(temp_data.n_cols,1);
        Wvector = Wvector / arma::norm(Wvector,2);
        
        // 正交向量处理
        if(Iter == 0){
            IterVector = arma::eye<arma::mat>(Wvector.n_elem,Wvector.n_elem);
        }
        else{
            // BeginEig = W_Record(:,end);
            // IMatrix = eye(length(BeginEig));
            // IterVector = IMatrix - BeginEig*BeginEig';
        }

        if(Wvector.n_rows != Sb.n_rows){
            Wvector = Wvector.t();
        }

        // 起始寻优向量
        W_Record.col(0) = Wvector;
        // 计算初始结果
        // 这里计算目标向量与数据相乘之后的比值
        J.col(0) = ObjFun(Sb,Sw,Wvector,IterVector);

        int times = 0;
        for(int i=0; i<100; ++i){
            // 计算提取，决定下一次的步进迭代方向
            double Gradient = GreFunc(Sb,Sw,Wvector,IterVector);

            arma::mat w_temp(Wvector.n_elem, gamma.n_elem, arma::fill::zeros);
            arma::mat ObjValue(1, gamma.n_elem, arma::fill::zeros);
            for(int gamma_i=0; gamma_i<gamma.n_elem; ++gamma_i){
                w_temp.col(gamma_i) = Wvector * (1 + gamma.at(gamma_i) * Gradient);
                // 将迭代向量结果，带入目标函数，求出目标函数值
                w_temp.col(gamma_i) /= arma::norm(w_temp.col(gamma_i));
                arma::mat w_temp_col = w_temp.col(gamma_i);
                double ObResult = ObjFun(Sb,Sw,w_temp_col,IterVector);
                ObjValue.col(gamma_i) = ObResult;
            }
            // printf("ObjValue.index_max()=%lld\n",ObjValue.index_max());
            Wvector =  w_temp.col (ObjValue.index_max());
            Wvector /= arma::norm(Wvector);
            W_Record.col(i+1) = Wvector;
            // Wvector.print("Wvector:");

            
            // 计算迭代结果
            J(i+1) = ObjFun(Sb,Sw,Wvector,IterVector);

            if(std::abs(J(i+1) - J(i)) < eta){
                break;
            }
            if(i > 2 && (J(i+1) < J(i))){
                // W_Record.col(i+1).zeros();
                break;
            }

            // 迭代次数++
            ++times;
        }
        // printf("times=%d\n",times);
        // printf("Iter=%d\n",Iter);

        // 不考虑PCA
        Eigvector.col(Iter) = W_Record.col(times);
    }

    // 消除重复映射向量
    // Eigvector = arma::unique(Eigvector);

    #ifdef DEBUG
    Eigvector.print("Eigvector:");
    #endif

    // 清除临时变量
    temp_data.clear();

    return 0;

}

double EmoRecBase::ObjFun(arma::mat &Sb, arma::mat &Sw, arma::mat &w, arma::mat &IterVector){
    return arma::norm(w.t() * IterVector * Sb, 1) / arma::norm(w.t() * IterVector *Sw ,1);
}

double EmoRecBase::GreFunc(arma::mat &Sb, arma::mat &Sw, arma::mat &w, arma::mat &IterVector){
    return arma::norm(IterVector * Sb, 1) / arma::norm(w.t() * IterVector * Sb, 1) 
        - arma::norm(IterVector * Sw, 1) / arma::norm(w.t() * IterVector * Sw, 1);
}



